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We numerically solve the inhomogeneous Zerilli-Moncrief and Regge- Wheeler equations in the time 
domain. We obtain the gravitational waveforms produced by a point-particle of mass ^ traveling 
around a Schwarzschild black hole of mass M on arbitrary bound and unbound orbits. Fluxes of 
energy and angular momentum at infinity and the event horizon are also calculated. Results for 
circular orbits, selected cases of eccentric orbits, and parabolic orbits are presented. The numerical 
results from the time-domain code indicate that, for all three types of orbital motion, black hole 
absorption contributes less than 1% of the total flux, so long as the orbital radius rp(t) satisfies 
rp(t) > 5A/ at all times. 
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Q I Tightly bound binary systems consisting of a compact object of a few solar masses and a supermassive black hole of 
10^ — 10^ Mq are very promising sources of gravitational waves for space-based detectors such as LISA [1]. There is 
now strong evidence that most galaxies harbour a 10^ — 10^ Mq supermassive black hole in their centre [2], and that 
, they are likely surrounded by a large population of solar- mass compact objects that reside in the galactic cusp [3]. 

The motion of objects in the galactic cusp is governed by the gravity of the supermassive black hole, but they are 
also constantly scattered due to the presence of multiple compact objects. For a given compact object, this process 
' occurs until it settles on a highly eccentric orbit that is tightly bound to the central black hole. On such an orbit, 
the object passes very close to the black hole at periastron and it emits a significant amount of gravitational waves. 
Capture occurs for those orbits that are sufficiently eccentric and have a sufficiently small periastron (on the order of 
M) [4]. In these cases, orbital evolution is driven by emission of gravitational waves, and the binary strongly radiates 
gravitational radiation, until the final plunge of the compact object into the central black hole. 
• The question is then to determine the rate at which solar-mass compact objects are captured by the central black 
hole and how quickly the orbits decay by emission of gravitational waves. Because capture occurs when the time to 



O i evolve due to emission of gravitational waves is much smaller than the time to evolve due to diffusion and scattering, 
^ determination of the type of orbits for which capture occurs and estimate of capture rates are sensitive to the strength 
5-H i of gravitational wave emission. Current estimates of orbital parameters for which capture occurs and associated 
capture rates are based on the quadrupole approximation for the emission of gravitational waves [5] . Although this 
is well justified for large periastron, it is not a good approximation for highly eccentric orbits with small periastron, 
' those of interest for gravitational-wave astronomy. 
rS \ In this paper, we consider a situation in which the compact object has already been captured by a spherically 
■ symmetric central black hole, and calculate the correct, general relativistic, rates at which the system loses energy 
and angular momentum to gravitational waves. We consider three types of orbits: circular, eccentric and parabolic 
orbits. These calculations will then be used to refine capture rates estimate, but this will be left for future work. 

At this level of approximation, the internal dynamics of the small compact object are irrelevant. We treat it as a 
point-particle and base our calculations on first-order perturbations of a Schwarzschild black hole; this is appropriate 
in view of the small mass ratio involved. The gravitational waveforms produced by the orbital motion are obtained 
by solving the even parity Zerilli-Moncrief [8, 11] (ZM) and the odd parity Regge- Wheeler [7] (RW) equations. We 
work with Schwarzschild coordinates, for which both wave equations take the form 
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where r* = r + 2Mlog(r/2M — 1) is the usual tortoise coordinate, and Vi{r) is a potential defined in Eq. (A4) for 
both modes. Explicit definitions for the ZM and the RW functions are given in Eqs. (A2) and (A3) of Appendix A. 
The source term Sim {f-, t) is of the form 

Sim{r,t) = G{r,t)6[r ~ Vpit)] + F(r,i)J'[r - Vpit)], (1.2) 

where a prime denotes an r-derivative, rp{t) denotes the radial position of the particle as a function of time, and 
G{r,t) and F{r,t) are known functions of r and t once the orbital motion of the particle is specified; they are given 
by Eq. (A6) for the Zerilli-Moncrief equation, and by Eq. (A7) for the Regge- Wheeler equation. 



2 



Instead of Fourier decomposing Eq. (1.1) and solving in the frequency domain, we choose to integrate them in 
the time domain. The numerical method we use was first developed by CO. Lousto and R.H. Price [6], and later 
corrected by K. Martel and E. Poisson to yield second-order convergence [12]; it is a finite-difference scheme, based on 
the null-cones of the Schwarzschild spacetime, which incorporates the source term without approximating d[r — rp{t)] 
and 6'[r — rp{t)]. This method is advantageous compared to Fourier decomposition because of the need, in the case of 
highly eccentric orbits, to sum over a very large number of frequencies in order to obtain accurate results [9, 15]. As 
an added bonus, the time-domain method provides the Zerilli-Moncrief and Regge- Wheeler functions everjrwhere in 
the spacetime. For each multipole moment, information about the fluxes of energy and angular momentum at infinity 
and through the event horizon is obtained by a single numerical integration. 

Astrophysical black holes are very likely to be rapidly rotating and the assumption of spherical symmetry for 
the central black hole is unrealistic. However, removing this assumption would require a substantial revision of our 
numerical method. The source term for the Schwarzschild perturbation equations can be treated exactly because, 
by removing the angular dependence, the problem is reduced to integrating a one-dimensional partial differential 
equation. A divergent source term of the form of Eq. (1.2) then leads to a simple jump in the field at the particle's 
position, and this can easily be handled by finite-difference methods. For a rotating black hole, one is faced with 
the task of solving the inhomogeneous Teukolsky equation [13]. It is well known that this equation is not separable 
in the time domain, because the eigenvalues of the angular functions are frequency dependent. Insisting on working 
in the time domain leaves a two-dimensional partial differential equation to integrate. Unfortunately, a (5-function 
source no longer leads to a simple jump at the position of the particle: the field is now (logarithmically) divergent at 
this location. Standard finite-difference methods are inadequate to deal with this type of behaviour and cannot be 
used. The problem can be circumvented by smearing the particle around its position (for example by using narrow 
Gaussian functions instead of (5- functions). This eliminates the divergence in the source term and, consequently, in the 
field; standard finite-difference methods can then be applied. Such an approach has been used to obtain gravitational 
waveforms produced by a particle on an equatorial circular orbit of the Kerr black hole [14], but the error introduced by 
smearing the particle is difficult to ascertain. By specializing to Schwarzschild, comparison with the present work will 
allow such a determination; this is another important justification for the work presented here. With this application 
in mind, we consider orbits with a wide range of eccentricities and semi-major axis, and do not necessarily restrict 
ourselves to highly eccentric orbits. 

The paper is organized as follows. In Sec. II we describe the orbital parametrization of bound and marginally- 
bound geodesies of the Schwarzschild spacetime. In Sec. Ill A we provide a relation between the Zerilli-Moncrief and 
Regge- Wheeler functions and the gravitational waveforms at infinity and near the event horizon; from these relations, 
the fiuxes of energy and angular momentum can be calculated. In Sec. IIIB, we provide a discussion of numerical 
issues that limit the accuracy with which we can determine the fluxes. In Sec. IIIC, Sec. HID, and Sec. HIE, we 
present our results for the gravitational waveforms and fluxes for circular, eccentric, and parabolic orbits, respectively. 
In Sec. IV, we summarize our findings. Appendix A contains a brief summary of first-order black hole perturbation 
theory, while Appendix B contains a detailed derivation of the flux formulae at infinity (Appendix B 3) and through 
the event horizon (Appendix B4). 



II. ORBITAL PARAMETRIZATION 



Following C. Cutler et al. [15], we introduce p, the semi-latus rectum, and e, the eccentricity, as orbital parameters. 
They are defined so that the periastron and apastron are at pM/{l + e) and pM/{l — e), respectively. In terms of 
these parameters, the energy and angular momentum per unit mass of a point particle are 



~2 ^ (p-2-2e)(p-2 + 2e) 



p{p — 3 — e^) 



p — 3 — e 



2 ■ 



(2.1) 



For e = the periastron and apastron coincide, and the orbit is circular. In the interval < e < 1, the motion occurs 
between two turning points, while for e = 1, the apastron is pushed back to infinity and the motion is parabolic^. In 
all cases, stable orbits exist only if p > 6 -|- 2e. 



^ In analogy with Newtonian mechanics, we use the term "parabolic" for marginally bound orbits; they have e = 1 (B = 1), but the 
trajectories traced out are not parabolae, except in the limit p S> 1. 
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FIG. 1: In the left panel, we display the trajectories in the Xp-yp plane for a geodesic with e = 1, and p = 8.001. For this 
choice of parameters, the particle orbits the black hole approximately four times before leaving the central region. In the right 
panel, we display a e = 0.9 and p = 7.8001 geodesic. When the particle reaches the periastron, it orbits the black hole on a 
quasi-circular orbit for approximately six cycles. In both cases, the exact number of cycles is given by Eq. (2.5). 



The position of the particle at time t is given by the coordinates {rp{t), fpit), 9p = 7r/2). Inspired by the solution 
to the two-body problem in Newtonian mechanics, the radial position of the particle is expressed as 

rpix) = , (2.2) 

1 + e cos X 

where x is a parameter along the orbit. This is well behaved at the turning points (x = 0, tt), which facilitates the 
numerical integration of the geodesic equations for the time and angular coordinates. In terms of Xi these are [15] 

-^t = Mp^ (p-2-2e)V^(p-2 + 2e)V^ 

dx {p — 2 — 2ecosx)(l + ecosx)^(p — 6 — 2ecosx)^/^ ' 



d 



dx^'' — 6 — 2e cosx)^^^ 



(2.4) 



The first of these equation can be numerically inverted to yield x{t)\ knowledge of rp(x) and <fp{x) is then equivalent 

to rp{t) and (pp{t). 

The geodesic equations given by Eqs. (2.3) and (2.4) are integrated using the Burlisch-Stoer method [10], and 
we chose the initial conditions as follows. The gravitational waveforms are extracted as functions of time at a 
location r*^^. We take the initial moment t = —to < to be the one at which the particle is at periastron (x = 0, 
fpi—to) = Mp/{\ + e), and (pp{—to) = 0). Wc set to equal to the light travel time between the periastron and the 
observation point. Thus, radiation emitted at the initial moment will reach the observer at t « 0. 

This parametrization of the geodesic is suitable for bound and unbound orbits of the Schwarzschild spacetime; for 
e < 1, the parameter x can take any real value, whereas for e > 1, it is confined to —ir/e < x 1^ n/e. In this paper, 
we consider circular orbits, selected cases of eccentric orbits, and parabolic orbits {e — 1), but the code is capable of 
producing gravitational waveforms for any value of e. For any p and e, the particle orbits the central black hole a 
number N = A(^p/(27r) of times before moving out of the central region. Integrating Eq. (2.4) over one radial period 
yields [15] 



N='-, r^~ K f J , (2.5) 

where K{m) — J^^"^ dx{l — m sin^ x)^^^'^ is the complete elliptic integral of the first kind. To visualize the trajectories, 
we introduce Xp{t) = rp{t) / {2M) cos{(pp{t)) and yp{t) = rp{t) / {2M) sm{(pp{t)) . In Fig. 1, we display trajectories in 
the Xp-Up plane for p = 7.8001 and e = 0.9 (left), and p = 8.001 and e = 1 (right). In both cases, the number 
of times the particle orbits the central black hole is large. This is because p is close to the critical value 6 + 2e at 
which N diverges. In these cases, gravitational-wave emission is dominated by the quasi-circular portion of the orbit 
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near pcriastron. The total energy emitted is then well approximated hy E = NEdrcuiar, where N is the (divergent) 
number of orbits, P the period of a circular orbit at Vp = Mp/{1 + e), and Edrcuiar is the energy emitted by a particle 
on such an orbit; a similar approximation holds for L. 



III. WAVEFORMS, ENERGY AND ANGULAR MOMENTUM RADIATED 



To numerically evolve Eq. (1.1) initial conditions must be provided for the gravitational perturbations. The manner 
in which the initial configuration of the gravitational field influences the subsequent evolution has been studied 
previously for radial geodesies [12]. For bound geodesies, the motion is quasi-periodic and waiting a sufRciently long 
time eliminates the contribution from the initial conditions, which simply propagates away. For marginally-bound 
geodesies, we chose the initial position of the particle to be very far from the pcriastron. Far away from the black hole, 
the velocity of the particle is small and it takes much longer for the particle to reach pcriastron than for the initial 
gravitational-wave content to escape from the system. At the point where the emission of radiation is strongest, 
there is no trace left of the initial configuration of the gravitational perturbations. This allows us to completely 
avoid problems related to the choice of initial data for both bound and marginally-bound geodesies. We chose zero 
initial conditions for the gravitational perturbations, acknowledging that this is inconsistent (creating the particle 
from nothing violates energy- momentum conservation), but recognizing that artifacts of this choice disappear in time. 
Fluxes may then be computed reliably after waiting a sufficiently long time. 



A. Far-zone fluxes and black hole absorption 



We first provide a short summary of the relations between the Zerilli-Moncrief and the Regge- Wheeler functions 
and the radiative portion of the metric perturbation at infinity and at the event horizon, as well as flux formulae used 
throughout the paper. Complete summaries of first-order black hole perturbation theory and flux calculations arc 
relegated to Appendices A and B, respectively. In the radiation zone, the two gravitational-wave polarizations are 
related to the Zerilli-Moncrief and the Regge- Wheeler functions by 



h+-ih^ = ^J2J f^fl [^ZM{t) - 2t y* dtVw(i')) -2Y'"^{e, v). (3.1) 



Irn 

Similarly, when r 2M, the two gravitational-wave polarizations are given by 



^++'^>' = m^ V ttM. (^^^^*^ (3.2) 

im V ^ ^ ' 

In these equations, gF'™ are spherical harmonics of spin- weight s [17]. 

From Isaacson's stress-energy tensor for gravitational waves [16], as well as Eqs. (3.1) and (3.2), we calculate the 
energy flux in each multipole moment to be 

^oo,eh _ Jei^I^IV-^Mp ,/ + meven 

i WII^I'^^M'r ,Z + modd, 

and the angular momentum flux to be 

rcx> eh _/ 1^1^ "^^^"^^^ + ,^ + "ieven 



W;^ V'w / d^Rw +C.C. , Z + m odd. 



where a " " " over a quantity denotes complex conjugation, c.c. is the complex conjugate, and the I and m indices are 
implicit on ipzM and ij^RW In Eqs. (3.3) and (3.4), E"^ and denote the fluxes across a surface r = const. oo, 
while Ef^ and Lf^ denote the fluxes through a surface r = const. 2M. The fluxes at inflnity are calculated using 
the Zerilli-Moncrief and Regge- Wheeler functions extracted at r* = r*jjg, where r*^^^ is large and positive, while for 
the horizon fltixes, they are extracted at r* = r*j^, where r*j^ is large and negative. Once Eim and Lim are known, the 
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total fliixes are obtained by summing over all modes: 

oo / 

^oo.eh ^ Y^^r'"^, E^'""" = E^'''' + 2^2^^''^, (3.5) 

1=2 m=l 

CSO / 



m=l 



there is no m = contribution to the angular momentum flux, and the factor of 2 in front of E^^ and i^'** comes 
from folding the m < contributions over to m > (see Appendix A). In a slow- motion, weak- field approximation 
the quadrupole moment dominates {1 = 2 and m = 2) and the total energy and angular momentum radiated over one 
orbital period are [19] 

647r/x2 , 73 2 , 37 _^/2 



LQip.e-l = ^„'(l + leAp-\ (3.8) 
The average energy and angular momentum radiated per unit time, defined by performing an orbital average, are 

-«) = f(^r^(-s»^-S'*). 



5 M pV2 V 8 , 

For the binaries considered in this paper, the slow-motion and weak-field approximations break down, and the fluxes 
must be computed using Eqs. (3.5) and (3.6). Numerically we cannot perform the infinite sums, and we truncate 
them at a finite value Zmax- In the next subsection, we explain the criteria used to choose /max and discuss the overall 
accuracy of the time-domain computation. 



B. Accurate determination of the fluxes: numerical issues 



In order to calculate the fluxes to a relative accuracy e (we use e = 0.01), we need to consider three sources of 
error: discretization of Eq. (1.1), effects of the finite size of our computational grid, as well as truncation of the sums 
in Eqs. (3.5) and (3.6). 

Firstly, discretization of Eq. (1.1) introduces numerical truncation errors. In a previous paper, we showed that the 
Lousto and Price algorithm can be corrected to yield second-order convergence [12], i.e. truncation errors scale as At"^, 
with At denoting the numerical stepsize for evolution. Throughout this work we generated gravitational waveforms 
by setting At = 0.1 (2M) in the numerical algorithm; this proved sufficient to determine the fluxes at infinity to the 
desired 1% accuracy. However, for a given stepsize the fluxes through the event horizon are never determined as 
accurately as the fluxes at infinity. The gravitational waves flowing through the event horizon are weaker than the 
ones escaping to infinity, and, because of this difference in scales, horizon fluxes are determined with an accuracy 
< 5%. But we will see below that horizon fluxes never amount to more than a few percents of the total fluxes. The 
lower accuracy with which black hole absorption is determined is then siifficient for our goal of 1% overall accuracy. 

Secondly, the expressions for the fluxes displayed in Eqs. (3.3) and (3.4) hold only asymptotically (r* ±oo). 
Numerically we are forced to extract the waveforms at finite r* values, and this introduces finite-size effects in our 
results. Numerical efficiency requires a small computational grid, but accuracy requires the waveforms to be extracted 
at a large and positive r*j,g and at a large and negative r*^. The flux formulae developed here are based on the stress- 
energy tensor for gravitational waves, as constructed by Isaacson [16]. The validity of the construction depends on 
X/TZ <C 1 being satisfied, where A is a wavelength of the radiation and TZ a typical radius of curvature. To calculate the 
fluxes far from the black hole, we extract the waveforms in an approximate radiation zone defined by A/robs *C 1, where 
A-i - (M/i?3)i/2 and Rp is a typical orbital radius. The radiation zone is then defined by R.p/rohsiRp/Mf/'^ < 1. 
For relativistic motion Rp ^ M and by imposing Rp/vohs < £, we make an error of order e in approximating the 
radiation zone. This is somewhat different from the criteria for the validity of Isaacson's stress-energy tensor, but 
since TZ~^ ~ (M/r^i^^)^/^, we have that X/TZ ~ {Rp/robs)"^^^ ^ e'^/^, and the use of the stress-energy tensor is justified. 
In practice we also imposed r*^^ > 750(2M). At the horizon, the situation is somewhat different. The typical radius 
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of curvature is 7^ ^ \/22M, but the radiation is blueshifted so that A ^ /oh(-Rp/M)^/^ 0, where /ch = 1 — 2M/rch- 
The requirement X/TL < 1 then translates to /eh(-Rp/(2M))3/2 ^ g_ ^Ve used -R;/|r*h| < e, as well as r^^ < -750(2M), 
which amply satisfies the above requirement. This yielded good results, but a better, more efficient choice would have 
been r^h = 2M[1 + {2A^/RpY^'^e]. With these choices of r*^^ and r*^, we are making an error of at most order e in 
determining the fluxes at infinity and through the event horizon, respectively. 

Finally, the last source of error limiting the accuracy of the determination of the fluxes arises from truncating the 
sums in Eqs. (3.5) and (3.6) at a finite value Zmax- The error is made small enough for our requirements by demanding 
that 

/ E?° L?° \ 

e = max -N^^ < 1% (3.10) 

be satisfied. Typically, Eq. (3.10) is satisfied with Imax given by (p/(l + e)) ^ ' < e, which is known to hold 
for circular orbits [18]. Note that because Ef^^^ and Lf^^^ are included in the sum, the error comes from neglecting 
terms starting at I = Imax + 1- In effect, the relative error made from neglecting these terms is much smaller than 
1%. In the sequel, we will return with empirical estimates of our numerical errors; these will confirm the preceding 
qualitative discussion. 



C. Circular orbits 



For circular orbits, e = and the radius of the orbit is Vp = pM . In Fig. 2 we display typical gravitational waveforms 
emitted by a particle traveling on a circular orbit. Both waveforms have the same pattern: The field oscillates with 
an angular frequency given by mO, where fl = M~^p~^/'^ is the orbital angular velocity and m is the multipole 
index. The left panel contains the dominant quadrupolar mode {1 = 2 and m = 2), while the right panel contains the 
dominant odd parity mode [1 = 2 and m = 1). 




t/2M t/2M 



FIG. 2: The dominant radiation modes for the Zerilli-Moncrief (left, I = 2 and m = 2) and Regge- Wheeler (right, I = 2 and 
m = 1) functions for a particle orbiting the black hole at = 12M. At early times, the waveforms are dominated by the initial 
data content. We calculate the energy and angular momentum fluxes after a time t/{2M) = 350.0 

The code outputs Egr and Lgr directly, but it proves convenient to express the fluxes in terms of and cl'- 
coefficients that remain close to 1 for all values of p. The total fluxes are calculated using Eqs. (3.5) and (3.6) and we 
express the numerically obtained results in the form 

E^nip) = ceEq{p,0), 

L^j,{p) = CL Lq{p,0), (3.11) 

where Eq{p, e) and Lq{p, e) are given by Eq. (3.9) above with e = 0. For circular orbits, we should find E = QL and 
therefore, ce = cl- 

Circular orbits have been studied extensively and we use them to quantitatively test the accuracy of the time- 
domain method. We perform a comparison of our results with the time-domain code (TD) with results obtained 
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FIG. 3: In the left panel, we display ce{FD), as well as ce{TD) and cl{TD), as functions of p. Both ce and cl slowly approach 
1 from below for large p. For small values of p, the coefficients approach 1.15 as p approaches 6. In the right panel, we display 
the residuals Re and Rl as defined in the text. Using the time-domain method, the fluxes are calculated accurately to 0.7% 
for p = 6.0001, and to 0.2% for large values of p. 

in the frequency domain (FD) by E. Poisson [22]. In the left panel of Fig. 3 we display ce{TD), cl{TD), and 
ce{FD). In the right panel, we display the residuals, Re = IQQ\ce{TD)-ce{FD)\/ce{FD) and Rl = 100|cl(T£))- 
ce{FD)\/ce{FD). In the interval 6 < p < 50, the time-domain code reproduces the frequency domain calculations 
to 0.7% or better, with the best agreement occurring for large values of p. 

In Table I we perform a mode by mode comparison between the two methods for p = 7.9456 and p = 46.062. For 
p = 7.9456, the fluxes from each multipole moment calculated with the time-domain code agree to 1% or better with 
the fluxes calculated in the frequency domain. A similar agreement is found for p = 46.062, with the exception of 
the I = 3 and m = 1, and I = 4 and m = 1 modes, for which the relative difference is 2.4% and 4.1%, respectively. 
This results from the huge difference in amplitude between these modes and the dominant mode. The stcpsizc used 
throughout this work is sufficient to obtain an overall relative accuracy of 1%, but it is not small enough to determine 
individual, small-amplitude modes to better than 2 ~ 5%. Although these modes could be resolved properly by using 
a smaller stepsize, it is not necessary for our goal of 1% overall accuracy; the contributions from these modes to the 
total fluxes are six and ten orders of magnitude smaller than the leading-order contributions, respectively. As such, 
they do not affect the overall accuracy of the computation; for p = 7.9456 and p = 46.062 the total fluxes calculated 
with the time-domain method agree with the frequency domain results to within 0.2% and 0.3%, respectively. 

Black hole absorption was calculated in a weak-field and slow-motion approximation for a particle in circular orbit 
by E. Poisson and M. Sasaki [20] and K. Alvi [21] who showed that it gives rise to a correction to the quadrupole 
formula: E'^^/Eq = = L'^^/Lq, where v = p~^l'^ = {M/rpY/"^ is the orbital velocity. The time-domain method 
allows black hole absorption to be calculated for arbitrary geodesies. In particular, for circular orbits our results show 
that even wlic;n c ~ 0.4 and the particle travels in a region of strong gravitational field, the amount of energy and 
angular momentum absorbed by the black hole is always a small correction to the total fluxes. For highly relativistic 
motion, this never grows large enough to contribute more than 0.4% of the total fluxes (see left panel of Fig. 4). For 
the purpose of calculating total fluxes with an overall accuracy of 1%, black hole absorption can safely be ignored. 
The right panel displays the ratio of horizon fluxes to the fluxes at infinity, normalized by {M/rpY, the weak- field and 
slow-motion approximation. As expected for circular motion, the normalized ratio for energy and angular momentum 
are equal to each other, and they approach 1 for large p. 

We estimate the accuracy with which black hole absorption can be determined using the time-domain method to 
be 5%. This estimate is based on the following argument. For small amplitude modes, the accuracy with which 
their contribution to the total fluxes can be determined is limited by errors originating from the discretization of 
Eq. (1.1) and the finite stepsize used in the numerical evolution. Based on the accuracy of the / = 4 and m = 1 
mode for p = 46.062 in Table I, the error is seen to be < 5% for modes whose contribution is ten orders of magnitude 
smaller than the dominant mode. For the range of p values considered in this paper, black hole absorption is at most 
seven orders of magnitude smaller than the dominant contribution. (This is evaluated using the (M/vpY relation at 
Tp = 50 Af, the value at which black hole absorption is least significant.) It is then safe to assume that fluxes through 
the event horizon are determined with an accuracy < 5%. (For values of p close to 6 -|- 2e, black hole absorption is 
more important and therefore more accurately determined.) 
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TABLE I: Energy and angular momentum fluxes for circular orbits, calculated using a time domain (TD) code, are compared 
with fluxes calculated by E. Poisson usin^ a frequency domain (FD) approach [22]. Here we chose p = 7.9456 and p = 46.062. 
The energy fluxes are in units of (M/jj,) , and the angular momentum fluxes are in units of M/fj!^. They are calculated at 



r^bs = 1500M and r^ts 


= 5200M for p = 


7.9456 and p = 


46.062, respectively. 








I m 


E°° (FD) 


(TD) 


rel. diff. 


(FD) 


L°° (TD) 


rel. diff. 


p = 7.9456 


2 1 


8.1633e-07 


8.1623e-07 


< 0.1% 


1.82836-05 


1.8270e-05 


0.1% 


2 


1.7063e-04 


1.7051e-04 


< 0.1% 


3.8215e-03 


3.8164e-03 


0.1% 


3 1 


2.1731e-09 


2.17416-09 


< 0.1% 


4.86706-08 


4.8684e-08 


< 0.1% 


2 


2.5199C-07 


2.51646-07 


0.1% 


5.64396-06 


5.6262e-06 


0.3% 


3 


2.5471e-05 


2.54326-05 


0.1% 


5.7048e-04 


5.6878e-04 


0.3% 


4 1 


8.3956e-13 


8.3507e-13 


0.2% 


1.8803e-ll 


1.8692e-ll 


0.6% 


2 


2.5091e-09 


2.49866-09 


0.4% 


5.61956-08 


5.5926e-08 


0.5% 


3 


5.7751e-08 


5.7464e-08 


0.5% 


1.2934e-06 


1.2933e-06 


< 0.1% 


4 


4.7256e-06 


4.7080e-06 


0.4% 


1.0584e-04 


1.0518e-04 


0.6% 


5 1 


1.2594e-15 


1.25446-15 


0.4% 


2.82066-14 


2.80908-14 


0.4% 


2 


2.7896c- 12 


2.75876-12 


1.1% 


6.24796-11 


6.1679e-ll 


1.3% 


3 


1.09336-09 


1.08306-09 


1.0% 


2.44866-08 


2.4227e-08 


1.1% 


4 


1.2324e-08 


1.2193e-08 


1.1% 


2.7603e-07 


2.7114e-07 


1.8% 


5 


9.4563e-07 


9.38356-07 


0.8% 


2.11796-05 


2.0933e-05 


1.2% 


Total 


2.0317e-04 


2.02736-04 


0.2% 


4.54466-03 


4.5399e-03 


0.1% 


p = 46.062 


2 1 


1.8490e-ll 


1.8713e-ll 


1.2% 


5.7804e-09 


5.8497e-09 


1.2% 


2 


2.8650e-08 


2.87286-08 


0.3% 


8.95666-06 


8.98098-06 


0.3% 


3 1 


7.5485e-14 


7.7275e-14 


2.4% 


2.35986-11 


2.4158e-ll 


2.4% 


2 


1.0926e-12 


1.0990e-12 


0.6% 


3.4157e-10 


3.4359e-10 


0.6% 


3 


8.0640e-10 


8.08356-10 


0.2% 


2.52106-07 


2.52708-07 


0.2% 


4 1 


9.9792e-19 


1.03906-18 


4.1% 


3.11916-16 


3.24808-16 


4.1% 


2 


1.6018e-14 


1.6171e-14 


1.0% 


5.00756-12 


5.05558-12 


1.0% 


3 


4.6603e-14 


4.6799e-14 


0.4% 


1.4569e-ll 


1.46318-11 


0.4% 


4 


2.7937e-ll 


2.79976-11 


0.2% 


8.73396-09 


8.7525C-09 


0.2% 


Total 


2.9505e-04 


2.95846-08 


0.3% 


9.22396-06 


9.24868-06 


0.3% 



D. Eccentric orbits 

For eccentric orbits, < e < 1, and the radial motion is bounded by the periastron rp|min = pM/{l + e) and the 
apastron rp\max = pM/{\ — e). In Fig. 5 and Fig 6 we display waveforms for two cases: p = 12 and e = 0.2, as well 
asp = 7.801 and e = 0.9. 

This type of orbital motion generates gravitational waveforms that are different in nature and in frequency content 
from circular orbits. Rather than being emitted uniformly along the orbit, the radiation is now emitted preferably at 
periastron. As the eccentricity increases the radiation is emitted in short bursts occurring near periastron. In these 
situations a time-domain approach is far more efficient than a frequency-domain approach. The reason is that in order 
to correctly calculate the waveforms in the frequency domain, a large number of individual frequencies (harmonics of 
the radial and azimuthal frequencies) are required, and summing over them can be hugely expensive. By contrast, a 
time-domain method handles all frequencies simultaneously. 

The fluxes are calculated over a number of wave cycles according to 

<E> = 1^ Edt, (3.12) 

where T is a few (> 3) radial periods; a similar expression holds for < L >. To obtain a quantitative idea of the 
relative accuracy of the time-domain method for eccentric orbits, we compute the fluxes for two points in the p-e 
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FIG. 4: We display the energy and angular momentum fluxes through the event horizon normalized by the fluxes in the radiation 
zone. Even for highly relativistic motion, the horizon fluxes contribute less than 0.4% of the total fluxes. For circular orbits, 
the theoretical prediction is that E"^ /E'^ = L"^ /L'^ . Numerically, this relation is only approximate, but nevertheless the two 
curves arc indistinguishable. The right panel displays these ratios normalized by (rj,/M)~*, the weak-field and slow-motion 
approximation. 





FIG. 5: The Zerilli-Moncrief (left, Z = 2, m = 2) and Regge- Wheeler (right, l = 2,m=l) functions for p = 12 and e = 0.2. As 
in the case of circular orbits, early times are dominated by the initial data content. 



piano and compare our calculations with Cutler et al. [15]: i) p = 7.50478 and e = 0.188917, and ii) p = 8.75455 and 
e — 0.764124. The results are displayed in Table II. For small eccentricities, e.g. case (i), the agreement is similar 
to the agreement achieved for circular orbits. For large eccentricities, e.g. case (ii), the agreement is ~ 2%. Because 
eccentric orbits in Schwarzschild are characterized by two incommensurate frequencies, the gravitational waveforms 
are quasi-periodic. By working in the frequency domain. Cutler et al. were able to formally average their fluxes 
over an infinite time. It is not, of course, possible to perform such an average in the time domain. Rather, for high 
eccentricities, the fluxes are averaged over a limited number of radial cycles (~ 3). This difference in averaging the 
fluxes is the most likely source of disagreement between time-domain and frequency- domain calculations for case (ii). 
For case (i), this is not as much of an issue, since the radial period is short enough to allow the time average to be 
performed over 10 cycles or more. 

Finally, we calculate the total energy and angular momentum emitted during one radial period as functions of p for 
e = 0.5. We express the total energy and angular momentum radiated to infinity, as calculated from the time domain 
code, as 



EGR{p,e) = CE[EQ{p,e) + {N -l)EQ{p/{l + e),Q)], 
LGR{p,e) = CL[LQ{p,e) + {N-l)LQ{p/{l + e),0)], 



(3.13) 
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FIG. 6: The Zerilli-Moncrief (left, l = 2,m = 2) and Regge- Wheeler (right, l = 2,m=l) functions for p = 7.801 and e = 0.9. 
As in the case of circular orbits, early times are dominated by the initial data content. The radiation occurs in short bursts 
when the particle approaches the periastron. This is typical of the zoom- whirl behaviour studied in [9] . 



TABLE II: Comparison of averaged fluxes for eccentric orbits with Cutler et al. for two points in the p-e plane [15]. The two 
cases presented are: i) p = 7.50478 and e = 0.188917, and ii) p = 8.75455 and e = 0.764124. 





case 




Cutler et al. 


time domain 


rel. diff. 


i) 


p = 7.50478 


<E> 


3.1680C-04 


3.1770e-04 


0.3% 




e = 0.188917 


<L> 


5.9656e-03 


5.9329e-03 


0.5% 


ii) 


p = 8.75455 


<E> 


2.1008e-04 


2.1484e-04 


2.3% 




e = 0.764124 


<L> 


2.7503e-03 


2.7932e-03 


1.6% 




FIG. 7: Coefficients ce and cl for the energy and angular momentum radiated as functions of p for e = 0.5 eccentric orbits. 
Near the last stable orbit (p = 7), ce approaches 1.24, while cl approaches 1.26. 



where we use Egr = P{p, e) < E >, P{p, e) is the radial period of the orbit obtained by integrating Eq. (2.3) over 
< X < 27r, < S > is given by Eq. (3.12), N = N{p, e) is given by Eq. (2.5) with e = 0.5, Eq{p, e) and Lq{p, e) are 
given by Eq. (3.9), and Ce and Cl are parameters that stay close to 1 for all values of p. In Fig. 7 we display Ce and cl 
as functions of p for e = 0.5. The coefficient Ce is close to 0.9 for large p and approaches 1.24 for p near 7. Similarly, 
the coefficient stays close to 0.95 for large p and approaches 1.26 for p near 7. The formulae above for the total 
energy and angular momentum radiated by a particle in eccentric orbit are justified by the fact that they have the 
correct limiting behaviour both for p large and for p — > 6 + 2e. For large p, the total energy and angular momentum 
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radiated by a particle on an eccentric orbit are well approximated by the quadrupole formulae given by Eq. (3.7). In 
this limit, iV — s- 1 and Eq. (3.13) produces the correct approximate energy and angular momentum radiated. When 
p — > 6 + 2e, the particle orbits the black hole for a number iV — 1 of quasi-circular orbits whose radius is equal to the 
periastron radius Vp = Mp/{1 + e). In this limit is large and the second term of Eq. (3.13) dominates the energy 
and angular momentum radiated. This term corresponds to the energy and angular momentum radiated after N such 
quasi-circular orbits 




FIG. 8: Black hole absorption for a particle in an e = 0.5 eccentric orbit. The absorption of both energy and angular momentum 
is negligible until the particle reaches p w 7.3, at which point it contributes approximately 1 ~ 2% of the total fluxes; these are 
eccentric orbits whose periastron is smaller than 4.9M. The right panel displays the same ratio normalized by (M/rp|min)*- 

The frequency of the radiation emitted by the orbiting particle increases as the periastron of the orbit becomes 
smaller. Since for a given eccentricity e, the periastron is proportional to p, the frequency of the radiation increases 
with decreasing p. Because the potential barrier around the Schwarzschild black hole is less opaque to high-frequency 
gravitational waves, we expect an increase in black hole absorption with a decrease in p. This is confirmed numerically 
for e = 0.5 and displayed in Fig. 8. For p < 7.3 {r,p ss 4.9M), the absorption of energy and angular momentum by 
the black hole contributes more than 1% of the total fluxes, while for p > 7.3 it contributes less than 1% and can 
be ignored when determining the total fluxes. In the right panel of the figure, we display black hole absorption for 
e = 0.5, normalized by (M/rp|min)'*, where rp|min = Mp/{1 + e) is the periastron distance. This is the correction 
expected from black hole absorption for a particle in circular orbit at rplmin- We use this normalization here because 
black hole absorption for generic orbits has not been calculated analytically. For large p, black hole absorption for 
e = 0.5 does not seem to converge toward the slow-motion and weak-field approximation for circular orbits. The 
normalized energy stays above 1, while the normalized angular momentum curve stays below 1. But because the 
relation dE = fldL used in deriving black hole absorption for circular orbit does not hold in general, there is no reason 
to believe that (M/rp|min)** should hold for generic orbits. Determining the differences in black hole absorption due 
to a finite eccentricity in a weak-field and slow-motion approximation would require a more detailed analysis than 
ours, since it is in this regime that our determination of black hole absorption is the least accurate. 

For radiation emitted by a particle whose orbital parameters are p = 6 + 2e and < e < 1, the argument relating 
black hole absorption to the orbital separation suggests that black hole absorption should be an increasing function 
of e along the line p = 6 + 2e (rplmin is a decreasing function of e along this line). It then comes as no surprise that 
numerical results displayed in Fig. 9 support this assertion (we used p = 6.001 4- 2e). Along this line, the radiation is 
emitted principally at periastron, where the orbit is quasi-circular. The relation dE = fldL, where fl = (M/r^ jmin)^''^ 
is the angular velocity of a particle on a circular orbit at rp |min, holds approximately and we find E^^/E"^ k. L^^/L°°. 

E. Parabolic orbits 

Particles on a parabolic trajectory have e = 1 (equivalently E = 1), and p specifies the value of the periastron: 
''pimin = Mp/2 with p > 8. For large values of p, the particle does not spend much time around rplmin, the position 
where the radiation is maximum; the waveforms have a simple structure around f = 0, time at which the radiation 
emitted at rplmin reaches an observer at r^j^g. This is displayed for even and odd modes in Fig. 10. In contrast, when 
p approaches its minimum value (pmin ^ 8), the particle circles the black hole for a number N of cycles. Because N 
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FIG. 9: Displayed are E" /E°° and V /L°° as functions of eccentricity along the curve p = 6.001 + 2e. Because of the decrease 
in periastron distance with increasing e, black hole absorption increases with e. A good approximation to these curves is given 

by E'VB"" = (1 + lAe^) (i;<='^/£;°°) |e=o. 




t/2M t/2M 



FIG. 10: Displayed are i/'zm (left, I = 2,m = 2) and V'RW (right, Z = 2,m = 1) as functions of time for e = 1 and p = 40. As 
in the case of circular orbits, early times are dominated by the initial data content. Total energy and angular momentum are 
calculated between -300 < t/2M < 300. 

diverges at p = 8, we get the zoom- whirl behaviour displayed in Fig. 1 [9]. The quasi-circular nature of the motion 
when Tp approaches rp|,nin results in a number of oscillations in the waveforms; these occur near t = Q for the observer 
at r*ijg, and are displayed in Fig. 11 for p = 8.001. 

In similarity with eccentric orbits, we express the numerically calculated energy and angular momentum radiated 

as 

EGR{p,e) = CE[EQ{p,e) + {N -\)EQ{p/{\ + e)M, 

Lgr{p, e) = CL [Lq{p, e) + (iV - 1)Lq(p/(1 + e), 0)] , (3.14) 

where N — N{p,e) is given by Eq. (2.5) with e = 1, EQ{e,p) and LQ{e,p) are given by Eq. (3.9), and Ce and Cl 
are again parameters that stay close to 1 for all p. For parabolic orbits, the total energy and angular momentum are 
computed using Eqs. (3.5) and (3.6) as 

E^^ = f^E°°dt, 
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FIG. 11: Displayed are V'zm (left, Z = 2, m = 2) and V'RW (right, Z = 2, m = 1) as functions of time for e = 1 and p = 8.001. 
Early times, where the choice of initial data dominates, are not displayed in order to make the t = region clearly visible. 
The energy and angular momentum fluxes are integrated between —300 < t/2M < 300 to obtain the total energy and angular 
momentum radiated. 




FIG. 12: Coefficients ce and cl for the total energy and angular momentum radiated as functions of p for a particle in parabolic 
orbit. Near p = 8, ce approaches 1.81 while cl approaches 1.84. 



for T large (we used T = 300(2M)). In Fig. 12, we display ce and cl ioi parabolic orbits. These quantities are close 
to 1 for large p, but increase above 1 as p approaches 8. Near this value of p, ce reaches 1.81, while cl approaches 
1.84. As for eccentric orbits, we find that for largo p the energy and angular momentum approach the values given by 
the quadrupolc approximation, but that for p close to 6 + 2e they arc better approximated by the energy and angular 
momentum radiated by a particle orbiting the black hole N times on a circular orbit of radius rp|min = Mp/(1 + e). 

The argument given previously for eccentric orbits holds true for parabolic orbits: when p > 8, black hole absorption 
is more important than for circular or eccentric orbit (sec Fig. 9). Our numerical results show that for p > 10, 
and L'^^ account for less than 1% of the total energy and angular momentum radiated, while for j? < 10 they can 
contribute as much as 5% of the total amounts (see Fig. 13). Hence, for p < 10, black hole absorption contributes a few 
pcrccnts of the total energy and angular momentum radiated and needs to be included in an accurate computation. 
Black hole absorption is not determined as accurately as the energy and angular momentum radiated to infinity, but 
the error we make in evaluating it is never large enough to spoil our goal of ~ 1% overall accuracy. 

For completeness, in Table III we display -B^/f'^ and L^j^^, the total energy and angular momentum radiated to 
infinity and through the event horizon, as returned by the time-domain code, for a wide range of p values. Based on 
the accuracy obtained for circular and eccentric orbits, we estimate that the total energy and angular momentum lost 
to gravitational waves are calculated to a relative accuracy of 1 ~ 2%. The actual accuracy is likely to be close to 
the accuracy achieved for circular orbits. The reason for this is quite simple. For parabolic orbits, there is no issue of 
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FIG. 13: Black hole absorption for a particle following a parabolic geodesic. The absorption of both energy and angular 
momentum is negligible until the particle reaches p w 10 or rplmin ~ 5Af. The right panel displays the same ratio normalized 
by (M/rp|min)'', where rp|min is the radii at periastron. Again, this factor is meaningful only for circular orbits, and is used 
only to illustrate the behaviour of black hole absorption as a function of p. 

performing a time- average, since the particle passes through periastron only once and we calculate the total energy 
for that motion. 



IV. CONCLUSION 

The time-domain method can produce waveforms and compute the associated fluxes of energy and angular mo- 
mentum to a relative accuracy of a few percents. For circular orbits, the method is extremely reliable and produces 

fluxes with an overall accuracy of 1% or better over the whole range of p values explored. For eccentric orbits, the 
comparison with Cutler et al. [15] is spoiled by the difficulty in performing a time average of the fluxes over a suffi- 
ciently long time. Because the disagreement arises from the differences in time averaging, the time-domain method 
is still capable of producing accurate waveforms for highly eccentric motion. We stress here that the limitation is 
in the computation of the time-averaged fluxes, not in obtaining the waveforms. On the other hand, for geodesies 
with small eccentricities there is no such limitation and the time-domain results are in better agreement with those 
calculated by Cutler et al. [15]. In all cases, the time-domain method is capable of determining the fluxes accurately 
to 1 ~ 2%. Similar accuracy is obtained for the total energy and angular momentum radiated by a particle traveling 
on a parabolic orbit. 

We also computed the absorption of energy and angular momentum by the black hole. For circular orbits with 
p > 6, this contribution can always be neglected, but not for orbits whose periastron is smaller than 5M. For such 
orbits, black hole absorption contributes more than 1% of the total fluxes and cannot be ignored. We showed that for 
e = 0.5 it can constitute a correction as large as 2% of the total fluxes; for parabolic orbits the contribution increases 
to 5%. 
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APPENDIX A: A SHORT REVIEW OF BLACK HOLE PERTURBATION THEORY 

Since the pioneer work of T. Regge and J. A. Wheeler [7] and F.J. Zerilli [8], perturbations of the Schwarzschild 
black hole have been studied extensively. Here we provide a short summary of the formalism, including the source 
terms appropriate for the Zerilli-Moncrief and Regge- Wheeler equations. 
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TABLE III: Total energy and angular momentum radiated by a particle orbiting a Schwarzschild black hole in a parabolic 
orbit. As usual, Eqji and Lg"^ denote the energy and angular momentum radiated to infinity, while Eq^j^ and L'qjj are the 
energy and angular momentum absorbed by the black hole. The black hole absorption contributes to less than 1% when p > 10: 
for parabolic orbits with pcriastron smaller than 5M, black hole absorption contributes a significant amount to the total energy 
and angular momentum radiated. 



p 


771 OO 

^GR 


ZTieh 


TOO 

J^GR 


r eh 
J^GR 


8.00001 


3.6703E+00 


1.8876E-01 


3.0133E+01 


1.5208E+00 


8.001 


2.2809E+00 


1.1260E-01 


1.9088E+01 


9.1166E-01 


8.201 


7.1130E-01 


2.6586E-02 


6.6010E+00 


2.2142E-01 


8.401 


5.1740E-01 


1.6534E-02 


5.0433E+00 


1.4244E-01 


8.601 


4.0970E-01 


1.1376E-02 


4.1665E+00 


1.0148E-01 


8.801 


3.3767E-01 


8.1988E-03 


3.5706E+00 


7.5175E-02 


9.0 


2.8419E-01 


6.0880E^03 


3.1196E+00 


5.7026E-02 


9.2 


2.4409E-01 


4.6044E-03 


2.7756E+00 


4.3891E-02 


9.4 


2.1228E-01 


3.5352E-03 


2.4973E+00 


3.4211E-02 


9.6 


1.8644E-01 


2.7473E-03 


2.2665E+00 


2.6951E-02 


9.8 


1.6506E-01 


2.1565E-03 


2.0718E+00 


2.1428E-02 


10.0 


1.4712E-01 


1.7072E-03 


1.9048E+00 


1.7176E-02 


10.5 


1.1292E-01 


9.8226E-04 


1.5752E+00 


1.0193E-02 


11.0 


8.8979E-02 


5.8626E-04 


1.3320E+00 


6.2819E-03 


11.5 


7.1545E-02 


3.6067E-04 


1.1455E+00 


3.9971E-03 


12.0 


5.8467E-02 


2.2778E-04 


9.9827E-01 


2.6152E-03 


12.5 


4.8424E-02 


1.4732E-04 


8.7936E-01 


1.7546E-03 


13.0 


4.0567E-02 


9.7321E-05 


7.8158E-01 


1.2036E-03 


13.5 


3.4326E-02 


6.5582E-05 


7.001 lE-01 


8.4286E-04 


14.0 


2.9303E-02 


4.4999E-05 


6.3136E-01 


6.0123E-04 


14.5 


2.5134E-02 


3.1409E-05 


5.7133E-01 


4.3635E-04 


15.0 


2.1774E-02 


2.2273&05 


5.2088E-01 


3.2171E-04 


16.0 


1.6636E-02 


1.1675E-05 


4.3867E-01 


1.8241E-04 


17.0 


1.2978E-02 


6.4484E-06 


3.7499E-01 


1.0869E-04 


18.0 


1.0303E-02 


3.7241E-06 


3.2454E-01 


6.7569E-05 


19.0 


8.3029E-03 


2.2354E-06 


2.8383E-01 


4.3565E-05 


20.0 


6.7794E-03 


1.3882E-06 


2.5047E-01 


2.8990E-05 


22.0 


4.6735E-03 


5.8355E-07 


1.9949E-01 


1.3895E-05 


24.0 


3.3426E-03 


2.6965E-07 


1.6280E-01 


7.2550E-06 


26.0 


2.4638E-03 


1.3447E-07 


1.3549E-01 


4.0539E-06 


28.0 


1.8620E-03 


7.1373E-08 


1.1459E-01 


2.3931E-06 


30.0 


1.4374E-03 


3.9949E-08 


9.8213E-02 


1.4781E-06 


32.0 


1.1298E-03 


2.3372E-08 


8.5141E-02 


9.4827E-07 


34.0 


9.0223E-04 


1.4201E-08 


7.4534E-02 


6.2832E-07 


36.0 


7.2956E-04 


8.9299E-09 


6.5744E-02 


4.2804E-07 


38.0 


5.9799E-04 


5.7886E-09 


5.8484E-02 


2.9873E-07 


40.0 


4.9549E-04 


3.8422E-09 


5.2369E-02 


2.1296E-07 


42.0 


4.1455E-04 


2.6163E-09 


4.7168E-02 


1.5469E-07 


44.0 


3.4987E-04 


1.8214E-09 


4.2706E-02 


1.1427E-07 


46.0 


2.9763E-04 


1.2887E-09 


3.8849E-02 


8.5683E-08 


48.0 


2.5501E-04 


9.3316E-10 


3.5491E-02 


6.5130E-08 


50.0 


2.1993E-04 


6.8343E-10 


3.2550E-02 


5.0118E-08 



The perturbations of the Schwarzschild spacetime are described by a linear perturbation tensor h^i, = g^i, — 
gSchwarzschiid^ where Is the metric of the perturbed spacetime and gSchwarzschiid ^j^g Schwarzschild solution. This 
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tensor is written as a multipole expansion using scalar, vectorial, and tensorial spherical harmonics. In Schwarzschild 

coordinates, we have 
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where capital roman indices run over the angular coordinates {6, (p), f = 1 — 2M/r, F'"* are the usual scalar harmonics, 
in term of which the vectorial and tensorial spherical harmonics are defined as 
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Here, VLab ~ {^i^vc?0), a bar denotes the covariant derivative compatible with ^ab, and eab is the Levi-Civita 
tensor on the unit two-sphere. The coefficients of the multipole expansions have implicit I and m indices, and there 
is an implicit summation over these indices. Defined this way, Hq, Hi, H2, qo, qi, K and G are even parity modes, 
while /lo, h\ and /12 are odd parity modes. 

By construction, the vectorial and tensorial spherical harmonics obey 5'^"™ = (— )™S''™*, where S is any spherical 
harmonic function. (This relation holds for scalar spherical harmonics, and since vectorial and tensorial spherical 
harmonics are obtained by the action of real operators on y'™, it also applies to these functions.) An important 
consequence of this relation is that for real metric perturbations, the multipole moments must satisfy M''"™ = 
(— )™M''"*, where M'"' is any one of i/o, Hi, H2, qo, qi, K, G, /iq, hi and It is then easily established that 
^im^*im ^ ]^i-m.j^*i,-m_ rpj^jg justifies folding the TO < terms over to m > in Eqs. (3.5) and (3.6). 

A gauge transformation can be used to eliminate four of the metric perturbations. In the Regge- Wheeler gauge, 
this freedom is used to set qo = qi = G = 0, and /i2 = 0. The two scalar fields 



A ■ 1 



/ 



d 



K+^{ H2-r—K 



RW 



--hi, 
r 



(A2) 



(A3) 



where A = (Z + 2)(Z — l)/2 and A = A + 3M/r, are the Zerilli-Moncrief and Regge- Wheeler functions, respectively. 
Their evolution is governed by Eq. (1.1) with the potentials 



VzM 



Vrw = 



f 



r2A2 



2A2 A + 1 + 



6M 



+ 



ISAf-^ 



A+ — 
r 



and the source terms 



SzM = 




or 



+ r{A-f)Q"' + rfQ^ 



rA 



A(A - ly + (4A - 9)Mr + 



2/ 
A 



0'' 



'RW 



r ) dr 



These are constructed from the perturbing stress-energy tensor T^'^ and we have defined 

IG-Trr^ 



g«6 = 8Tr J T"' 



327rr4 



{l-l)l{l + l){l + 2) 



rj^AByhn* 



(A4) 
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and 



pa 



1 1 

WTT) 



rpaAj^lm* 



dfl, 



P = 



167rr 



lower-case roman indices run over t and r, the integration is over the unit two-sphere with dVl = sinO d9d(p, and T"*", 
T""^ and T^^ are components of T'*''. For a particle traveling on a geodesic with proper time r, coordinates z^{t), 
and four- velocity u^{t), the stress-energy tensor is 



(A5) 



where (5^ [a;" — 2"] is a four-dimensional Dirac functional. 

These expressions for the source terms can be used, combined with the geodesic equations of the Schwarzschild 
spacetime, to calculate explicit expressions for the factors G{t,r) and F{t,r) appearing in Eq. (1.2). For even parity 
modes (ZM), we get 



G(r, t) = a Y*{t) + b z;{t) + c u;^{t) + d v;^{t) 



(A6) 



where V'^ = f(l + P/r"^^ , 




A + 1- 



3M 



and 



IGtt L f 



Stt L2 f3 



Finally, the source terms for odd parity modes (RW) are 



X + 3 



-327r 



ii-2y. L^ p 

(/ + 2)! E r3 



G{r,t) 
F{r,t) 



-a- 



4 
r 



1 



3M 
r 



(A7) 



where 



A 



Note that G{r,t) and F{r,t) contain scalar, vectorial and tensorial harmonic functions evaluated at the angular 
position of the particle <fp{t). For example, for even modes, some terms in G and F are proportional to y''"(7r/2, (pp{t)). 
Because the orbital motion takes place in the equatorial plane, each spherical harmonic function is evaluated at 
9p = 7r/2. A useful consequence of this is that the source term for the Zerilli-Moncrief function vanishes when I + m 
is odd, while the source term for the Regge- Wheeler function vanishes when Z -|- m is even. This was used in Eq. (3.3) 
and Eq (3.4). 

Once V'ZM and V'rw are found by solving Eq. (1.1), with the source terms of Eqs. (A6) and (A7), the perturbation 
tensor can be reconstructed. In the Regge- Wheeler gauge, we have 



K = /_^zM+>l(r)V'ZM-^^;^ 



H. 



2 = -r 



A + 1 



ZM 



K 



4k, 

or 



Hi 

Ho 



d 
r — 
dt 

H2 - 



d_ 

dr 



^ZM + B{r)ip2,M 



x + 1 



r/_9 

Tdt 



0* 



(A8) 
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for even parity modes, where we have defined A{r) = [A(A + 1) + 3M/r (A + 2M/r)]/(rA) and B{r) = [A(l - 3M/r) - 
3M^/r'^]/(r/A). For odd parity modes, the reeonstructcd metric perturbations are 



ho = -fj dt' 



hi = -r/-VRW- (A9) 

APPENDIX B: RADIATION ZONE FLUXES AND BLACK HOLE ABSORPTION 

The fluxes of energy and angular momentum can be obtained from Isaacson's stress-energy tensor for gravitational 
waves: 

where (. . .) denotes an average over a region of spacetime large compared with the wavelength of the radiation. 
Typically, T^^ can be defined when the wavelength of the radiation. A, is small compared to a typical radius of 
curvature TZ. By definition, A <C 7^ in the radiation zone and the stress-energy tensor for gravitational waves can be 
defined there. There is a second region where the condition X ^ TZ is satisfied: A stationary observer near r = 2M 
sees TZ ~ 2M, but the radiation is strongly blueshifted and A ^ 0; that this is the case is clear from the divergence 
in Eq. (B30) below. 

Because the Schwarzschild black hole is static and axially symmetric, it possesses two Killing vectors, and 
that can be used, in conjunction with T^,y, to obtain expressions for the fiuxes of energy and angular momentum 
through a surface S at constant r. We have that (dropping the "GW on the stress-energy tensor) 



dE = 
dL 



-J^nit)CdE^ (B2) 
/ T'i (^)r dS^, (B3) 



where dS^ is an outward oriented surface element on S. In Schwarzschild coordinates, these expressions reduce to 

E = -er'^f j dn Ttr , 

L = er^f j dil Tr^ , (B4) 

where an overdot indicates differentiation with respect to t, and e is 1 when calculating the fluxes in the radiation 
zone and -1 when calculating black hole absorption. Because the event horizon is a null surface, it is conceptually 
better to use dE/dv and dL/dv in place of dE/dt and dL/dt (with v = t + r*). Similarly, because radiation travels 
to i7+, which is another null surface, it is conceptually better to use dE/du and dL/du instead of dE/dt and dL/dt 
for outgoing fiuxes (with u = t — r*). However, because we numerically extract the gravitational waveforms at finite 
values of r*, without ever reaching the event horizon or ^7+, the use of dE/dt and dL/dt is a better representation of 
our numerical procedure. 

We provided a summary of the perturbation formalism in Appendix A. This included explicit formulae to reconstruct 

the metric perturbations in the Rcgge- Wheeler gauge. To construct Ttr and Tr^, it proves convenient to extract the 
radiative part of the perturbation tensor. To this end, we introduce the null tetrad l'^, n^, and m'^: and ny are null 
vectors tangent to outgoing and ingoing rays, respectively, while m'* is a complex null vector (with complex conjugate 
m'') on the two-sphere. They satisfy the relations 

= = n'^n^, m'^m^ = = fn^fhij,, l^n./^ = -1 = -m^fhu. (B5) 

In Schwarzschild coordinates, their components are 

= (-!,/-!, 0,0), n^ = -i(/, 1,0,0), m^ = ^r(O,O,l,zsin0). (B6) 

To find the perturbation tensor in a radiation gauge, we seek a gauge vector that generates the transformation 
between the Regge- Wheeler gauge and the radiation gauge. Under such an infinitesimal coordinate transformation. 
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the perturbation field transforms as = hj^^ — ■ where RG stands for radiation gauge and RW for Regge- 
Wheeler gauge. Using the spherical harmonic functions introduced previously, the vector can be expressed as a 
multipole expansion [7]: 

^^°'''') = {o,oyKX'r), (B7) 

where at, a^, /3 and k arc freely specifiable gauge functions. Combining this with the multipole expansion for the 
perturbation tensor in the Regge- Wheeler gauge, we obtain the transformation law for the each perturbation mode: 

KRG = i^RW _ (V_^^ _ ^ ^ 
for the transformation of even parity modes, and 

for odd modes; an overdot designates a t-derivative and a prime an r-derivative. In the next two sections, we develop 
solutions appropriate to outgoing and ingoing radiation gauges, respectively. 



1. Outgoing radiation gauge 

In an outgoing radiation gauge h'^^'^, the perturbation tensor satisfies [23] 

h';^^^n''m'' = = h^^^n^nf, 
h^l^^nn" = = h'^^'^'m^ffir (BIO) 

The first two conditions indicate the h^^}'^ is transverse to outgoing null rays, while the last condition indicates that 
it is traceless. Eqs. (BIO) involve five conditions, one too many for the specification of a gauge. But this system is 
not overdetermined: once four of these equations are enforced, the fifth is found to be satisfied automatically in the 
radiation zone. 

From Eqs. (B8) and (B9), the gauge conditions can be expressed in terms of multipole moments. To leading-order 

in r~^, we get 

4(d,-d.) = Hr + Hr-2Hr', 

at-ar + 2r^l3 = 0, 

at + otr = 0, 

-ar = -K^"^, (Bll) 
r 
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for even parity modes, and 

2r2/i = hf^-h^"^, (B12) 

for odd parity modes. We used d/dt = —d/dr + 0{r~^), appropriate in the radiation zone, to eliminate r-derivatives 
in favor of t-derivatives. 

The right-hand side of these equations contains the perturbation tensor in the Regge- Wheeler gauge. From Eq. (A8), 
we find that, in the radiation zone, even parity modes have the asymptotic form 

li-RW ^ ,/, rrRW ^ rj-RW ^ rj-RW ^ „,/: 

K pa —tpZM, il2 ~ -"0 ~ ~ rtpzM, 

while, from Eq. (A9), we find that odd parity modes are given asymptotically by 

ho « -hi w rV'Rw- 

Solutions to the gauge transformation are then easy to find. For even modes, the last of Eq. (Bll) yields ar = 

— (l/2)r'0zM- The other gauge functions are easily obtained from the remaining equations of Eq. (Bll): at = —ar 
and /? = — ■0zM/(2r). For odd modes, direct integration of Eq. (B12) yields k = 1/r dt' ip^^ [t' ) . Going back to 
Eqs. (B8) and (B9), we can reconstruct the perturbation tensor in the radiation zone. Solving these for G'^^'^ and 
h^^^, we get 

/^Si'' = r (^VzM(i)Fi'g - 2 y* dtVRw(t')^AB) + 0{l). (B13) 

In this expression, the Zerilli-Moncrief and Regge- Wheeler functions have I and m multipole indices and there is an 
implicit summation over them. Dropping the 0{1) term, we refer to h^^^ as the radiative part of the perturbation 
tensor in the radiation zone, because it contains all of the information about energy and angular momentum carried 
to infinity. This last expression can be re- written in terms of the two gravitational- wave polarizations, /i+ = h^p''^ /r^, 
and /ix = h'^^^ /{r'^sm.e). The result is Eq. (3.1). 



2. Ingoing radiation gauge 



To obtain the radiative part of the gravitational field in the vicinity of the event horizon, we impose an ingoing 
radiation gauge. We seek a solution to leading order in / ^ 0, the expansion parameter near the horizon. 

The ingoing radiation gauge can be obtained from Eq. (BIO), by making the replacement ^ of the tetrad 
vectors. The same comment about the number of gauge conditions can be made here: only four gauge conditions 
need to be imposed, and the fifth condition is then satisfied automatically. 

From Eq. (B8), Eq. (B9), and the ingoing radiation gauge conditions, we obtain 

4(dt + "r) - ^ {at + ar) = 2/ {H^ + H^^) , 

at - ar = 0, 
8M'^P + at + ar = 0, 

j^ar-p = K^'^, (B14) 

for the gauge transformation of even parity modes, and 

r^k = r{h^^ + fhf^), (B15) 

for odd parity modes: we substituted ar f~^ar in Eq. (B8), we used d/dr* = d/dt + 0{f) to eliminate derivatives 
with respect to r*, and an overdot denotes a time derivative. 

The asymptotic form of the metric perturbations in the Regge- Wheeler gauge is obtained from Eqs. (AS) and (A9). 
We get 

K^-^ « ^zM + ^^ZM, nr = = Hr = r' (smvJzm - ^^zm) , 

for even parity modes, and 

hr=fhr=-2Mi;nw, 
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for odd parity ones. These can be inserted baek into Eqs. (B14) and (B15). The solution to the gauge transformation 
then proceeds as foUow. For even parity modes, the second and third equations yield = at and AM'^fi = —at- 
These can be substituted into the first and the time derivative of the fourth of Eq. (B14), to yield a system of equation 
for at'. 



at 



A + 1 
2M 



at = Mil)', 



A + 1 



ZM 



V'; 



ZM- 



(B16) 
(B17) 



Eliminating the time derivative by subtraction, we find at = MipzM- For odd parity modes, integration of Eq. (B15), 
combined with the asymptotic form of and /if^, yields k = — 1/(2M) J dt'tpw/vit')- 

Prom these, and Eqs. (B8) and (B9), we obtain G^^^ and h^^'^ , which are used to reconstruct the gravitational 
perturbation tensor: 



2M 



^zmVX^ + 2 



dtVRw(i')<B 



+ o{f). 



(B18) 



Again, the Zerilli-Moncrief and Regge- Wheeler functions have / and m multipole indices, and there is an implicit 

summation over them. The components h'^j^ (without the 0(f) correction) contains all the information about the 
energy and angular momentum absorbed by the black hole. It is then meaningful to refer to these components as the 
radiative part of the perturbation tensor in the vicinity of the event horizon. In analogy with the far zone definitions, 
the two gravitational-wave polarizations are defined as h+ = /ig^'^/4M^ and hx = h^^^/ (4M^ sin 9). They are given 
in Eq. (3.2). 



3. Radiation zone fluxes 

In terms of its tetrad components, the perturbation tensor in the outgoing radiation gauge is 

+ hrhfhiri^m^ + h.mm'rniji.fhu-, (B19) 
where hw = h^i^v^v'^, for any vector belonging to the tetrad. To leading-order in r~^, the tetrad components are 

hii ~ 0(r-'), 
him ~ 0{r-^), 

hfnrh = ^h%'^ fh^ffi'^ , 

hmm = ^h%°m^m^, (B20) 

where the vector = r {d0^/dx'^) m'', and h'^§'^ is given in Eq. (B13). 

Calculating the covariant derivative of h^^*^ and substituting the result in Eq. (Bl), we get 



647r 



+ C.C. , (B21) 



where 



Paffn = /lH["a«/3];/i + 2/lJm[n(am^)];^ + 2/l(^[n(am^)];^ 

+ hrnm[fnafhi}];^i + hfnm[marn0]-fj,. (B22) 

The first term appearing in Eq. (B21) can be calculated exactly. It is 

V ^iiHajiv ~ hmm,nh^fj^ ii + hfnin^^h^^ ^,. (B23) 
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Evaluating the remaining terms requires more effort. They involve products of tetrad vectors and their covariant 
derivatives. It is easy to show that the only non- vanishing t, r, and ip components are m^riaj^ = — (•\/2/4)z/sin0, 
m^maitp = —icos6, and m".^na;r = \/2M/(4r^)isin^. Using these and Eq. (B20), we find that 



1 

"I" '^^{hmmh ffiff^ -\- hmfhh^^) COS 6 + C.C., 



(B24) 



(B25) 



where we replaced r-derivatives with f-derivatives, and neglected terms of order 0{r~^) and higher. 

These expressions for the stress-energy tensor can be used to calculate the fluxes in the radiation zone. Inserting 
Eq. (B24) into the first of Eq. (B4), where we set e = 1 and / « 1, yields 



dE r'' f ■ ■ ■ ■ 1 f ■ ■ 



dt 32'7r ^ 

Im I'm' 

{l + 2)\ 



AB 



—y 

647r^ {l-2)\ 

Ira 



l^i-ZMl^ + 4|^RW|^ 



(B26) 



where in the first line we use ^^''^^^^hAsh^jj = r^{hmmh'!^fh + ^mmh'^^), the second line follows from Eq. (B13), 
and the third line follows from evaluating the angular integral with the aid of 



1 {I + 2)! 

2 (/ - 2)! 



The angular momentum fiux calculation follows similar steps. Inserting Eq. (B25) into the second of Eq. (B4), we 



get 



dL 
dt 



647r 



/ 



dQ 



C.C. 



(B27) 



The last term involves the product of cos^ with a term of the form sin 9 S^"^ (6) S*,^, (6), where S^'^{9) is a spherical 
harmonic function. Under the interchange 9 ^ tt — 9, wc have that cos0 —cos9 and sin9S''™{9)S*,„^,{9) 
siu9S^"^{9)Si,^,{9). The overall term is therefore odd in 9 with respect to ■k/2, and integration over < 9 < tt yields 
zero contribution to the angular momentum flux. We are then left with 



dL 
'dt 



647r 



647r 



/ 



{hrnm^rh 



^Jdn {h^^h*AB + C.C.) 



im 



dn 



Im I'm' 

im {l + 2)\ 



+ C.C. 



Im 



12877 (/ - 2)! 



1pZM V'ZM + 4V'RW 



/* dt'rTmit') 

J — oo 



+ C.C. 



(B28) 



where wc use hrnm,(p = —I'fnhmm in the first line, and the remaining steps shadow the ones for the energy fiux 
calculation. 
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4. Black hole absorption 

The calculation of the black hole absorption is similar to the calculation of the far-zone fluxes. Here we are looking 

to isolate the divergent piece of T^,y, since it is this part that corresponds to the blucshiftcd gravitational waves: The 
expansion parameter is /, and we are looking for the 0{f~^) portion of the tr and r(p components of T^,^. We neglect 
terms of order 0{1). 

The material developed in Sec. B3 can be used here simply by replacing Z** <-> n^. In an ingoing radiation gauge, 
the non-trivial tetrad components of the perturbation tensor are 

hnn ~ 

hum ~ Oif), 



h-- - - A -B 

1 



hmm = j^h'^^ m^m'' , (B29) 

and was introduced previously. 

With the replacement ^ Z", the steps we follow are almost exactly the same as those of the radiation zone 
calculations. The stress-energy tensor is written as in Eq. (B21), with 770/3 and Papv changed to reflect the exchange 
of tetrad vectors. It is not difiicult to show that the non-vanishing components of the contracted derivatives of the 
tetrad vectors are now m"la;ip = {V^/2)ism9, fh"ma;ip = —icosO, and m°'.^la;r = V^M / {2r^)if~^ s\n6 . 

These, combined with Eq. (B29), reveal that the relevant components of the stress-energy tensor are given by 

_ /-I . 

^tr — ~ '^l^ij^mmhffijfi -\- hyfirhhjfij^ (B30) 

/-I . ^ • * 

Trip = ~'^^{hmmhff^,j^ + hjfimh^^^^) 

+ -^{hmmh^^ + hfnfhh;^^) COS e + c.c, (B31) 

where we replaced r*-dcrivatives with t-derivatives and neglected term of order 0(1). Note that this is exactly of the 
same form (apart from a factor of f~^) as that obtained for Ttr and T^^ in the far zone. 

To calculate the fluxes, we insert these expressions into Eq. (B4), where we also set e = —1. The divergence in the 
stress-energy tensor is canceled by the factor of / appearing in Eq. (B4). The remaining calculations are identical 
with those of the far-zone, with given by Eq. (B18). The energy flux is 



dE _ ^ I {1 + 2)1 
Itt ~ f^64^(/-2)! 

Im ^ ' 



|V'ZMr + 4|V'Rw|^ 



while the angular momentum flux is 



dL _ -s-^ rni {1 + 2)! 
'dt ~ 12877 {I - 2)! 



(B32) 



c.c. . (B33) 
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